Groundwater depletion in California’s Central Valley accelerates during megadrought

Groundwater provides nearly half of irrigation water supply, and it enables resilience during drought, but in many regions of the world, it remains poorly, if at all managed. In heavily agricultural regions like California’s Central Valley, where groundwater management is being slowly implemented over a 27-year period that began in 2015, groundwater provides two–thirds or more of irrigation water during drought, which has led to falling water tables, drying wells, subsiding land, and its long-term disappearance. Here we use nearly two decades of observations from NASA’s GRACE satellite missions and show that the rate of groundwater depletion in the Central Valley has been accelerating since 2003 (1.86 km3/yr, 1961–2021; 2.41 km3/yr, 2003–2021; 8.58 km3/yr, 2019–2021), a period of megadrought in southwestern North America. Results suggest the need for expedited implementation of groundwater management in the Central Valley to ensure its availability during the increasingly intense droughts of the future.


S1. Water fluxes for water balance analysis
The water flux components including precipitation (P), evapotranspiration (ET), and streamflow discharge (Q) from PRISM, PT-JPL, and USGS gauges were integrated to monthly totals and expressed as basin-averaged depths in millimeters (mm) for water balance analysis. The PRISM and PT-JPL datasets originally record P and ET as monthly cumulative values in mm, so their time series data were extracted by averaging values over the domain as shown in Fig 1. The Q values from USGS gauges were recorded in cubit-feet per second on a daily basis (Qcfs). To convert Qcfs to monthly values in mm (Qmm), the equation (S1) is used: where, d represents the d th day of the month and n is total days of the month. q1 is a conversion factor, 2.83168x10 -11 , for cubit-feet to km 3 , and A is the area of the study region, equivalent to 153,659 km 2 . The P, ET, and Q values are presented in Fig. 2(B) and used to calculate the water balance using equation (1) as shown in the Fig. 2(C) in the main text.

S2. Surface water components for groundwater storage change estimation
The monthly soil moisture data from NLDAS is recorded in kg/m 2 . Using conventional water density, i.e. 1000 kg/m 2 , the SM value of NLDAS is equivalent to Equivalent Water Height in millimeters (EWH, mm).
The surface water storage from in situ gauges of dams and reservoirs are recorded in units of acre-feet (SWAF). In the study, surface water storage (SW) is converted to EWH in millimeters (mm) using the equation below: where, q2 is a conversion factor, 1.23346x10 -6 , for acre-feet to km 3 , and A is the area of the study region, equivalent to 153,659 km 2 .
The snow water equivalent from SNODAS is recorded as EWH in meters and can be directly converted into EWH in mm. The three surface water components were extracted for the study region and the three sub-basins, and subtracting their historical means

S3. Uncertainty of groundwater recharge and depletion rates
The groundwater recharge or depletion rates (b) were calculated using linear regression with the linear equation as below: where y is the groundwater anomaly, x is the time step in years, and c is the intercept of the linear equation. We quantified the uncertainty of groundwater change rates (∆ ) by using student's t test.
where, t is the t score from statistic t table which depends on sample numbers and the definition of the confidence interval (95% in the study), and SEb is the standard error of the coefficient estimate.
where, 6 and ̂ are the slope (estimated groundwater change rate) and intercept from the regression, respectively; 6 and 6 are the groundwater anomalies at time x of sample i during the regression period; n is the total sample number for the regression, and DR is the degree of freedom, equivalent to n-2 for linear regression.

S4. Water table depth from wells
The in situ well measurements used in the study are identical to those used in Kim et al., (1), which were compiled and processed using groundwater monitoring networks managed by California's Department of Water Resources (DWR) and the U.S. Geological Survey (USGS). The archived raw data are inconsistent in well identification system, format, recording period and temporal resolution, and quality etc., which presents a challenge to integrate and merge the datasets for regional domain monitoring. Kim et al., (2021) (1) therefore developed a series processing scheme for data integration into a large domain. First, a filtering process was conducted based on Quality Control criteria provided from each management agency to remove erroneous measurements of each well. Then, these qualified measurements from all wells were merged and assigned a new identifier.
Because coordinates of wells from raw data do not have sufficient resolutions, some wells share the same or very close coordinates but with inconsistent values and nomenclature. To overcome this issue and only select dependable data, Kim et al., 2021 (1) applied a gridding and scoring system to only select high quality, unique data within each grid cell. The Central Valley domain was gridded into 1 km cells. All wells fell within an individual grid cell were scored based on 1. total number of measurements, 2. time coverage during the study period, and 3. mean measurement number per month. These criteria were normalized and weighted evenly to all wells. Measurements from the highest ranking well was then selected to represent the water table depth for that grid cell. Figure S2 shows the well stations distributed in the study region. The monthly water table depth of Central Valley and each basin domain were calculated by averaging the selected data across the domain, as shown in Fig. S3. Note that the well measurements from Sacramento basin may dominate water table depth variations for Central Valley due to about 4-6 times more samples as compared to the other two basins. Sample numbers for the entire Central Valley decrease by about 70% after 2019, which may result in an inconsistency in the time series of data, and may cause the discrepancy between well measurements and GRACE/FO estimates as shown in Fig. 3(B) in the main text.